home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / gaussfit.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  191 lines
  1. ; $Id: gaussfit.pro,v 1.6 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1982-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;
  6.  
  7. PRO    GAUSS_FUNCT,X,A,F,PDER
  8. ; NAME:
  9. ;    GAUSS_FUNCT
  10. ;
  11. ; PURPOSE:
  12. ;    EVALUATE THE SUM OF A GAUSSIAN AND A 2ND ORDER POLYNOMIAL
  13. ;    AND OPTIONALLY RETURN THE VALUE OF IT'S PARTIAL DERIVATIVES.
  14. ;    NORMALLY, THIS FUNCTION IS USED BY CURVEFIT TO FIT THE
  15. ;    SUM OF A LINE AND A VARYING BACKGROUND TO ACTUAL DATA.
  16. ;
  17. ; CATEGORY:
  18. ;    E2 - CURVE AND SURFACE FITTING.
  19. ; CALLING SEQUENCE:
  20. ;    FUNCT,X,A,F,PDER
  21. ; INPUTS:
  22. ;    X = VALUES OF INDEPENDENT VARIABLE.
  23. ;    A = PARAMETERS OF EQUATION DESCRIBED BELOW.
  24. ; OUTPUTS:
  25. ;    F = VALUE OF FUNCTION AT EACH X(I).
  26. ;
  27. ; OPTIONAL OUTPUT PARAMETERS:
  28. ;    PDER = (N_ELEMENTS(X),6) ARRAY CONTAINING THE
  29. ;        PARTIAL DERIVATIVES.  P(I,J) = DERIVATIVE
  30. ;        AT ITH POINT W/RESPECT TO JTH PARAMETER.
  31. ; COMMON BLOCKS:
  32. ;    NONE.
  33. ; SIDE EFFECTS:
  34. ;    NONE.
  35. ; RESTRICTIONS:
  36. ;    NONE.
  37. ; PROCEDURE:
  38. ;    F = A(0)*EXP(-Z^2/2) + A(3) + A(4)*X + A(5)*X^2
  39. ;    Z = (X-A(1))/A(2)
  40. ;    Elements beyond A(2) are optional.
  41. ; MODIFICATION HISTORY:
  42. ;    WRITTEN, DMS, RSI, SEPT, 1982.
  43. ;    Modified, DMS, Oct 1990.  Avoids divide by 0 if A(2) is 0.
  44. ;    Added to Gauss_fit, when the variable function name to
  45. ;        Curve_fit was implemented.  DMS, Nov, 1990.
  46. ;
  47.     n = n_elements(a)
  48.     ON_ERROR,2                      ;Return to caller if an error occurs
  49.     if a[2] ne 0.0 then begin
  50.         Z = (X-A[1])/A[2]     ;GET Z
  51.         EZ = EXP(-Z^2/2.)    ;GAUSSIAN PART
  52.     endif else begin
  53.         z = 100.
  54.         ez = 0.0
  55.     endelse
  56.  
  57.     case n of
  58. 3:     F = A[0]*EZ
  59. 4:     F = A[0]*EZ + A[3]
  60. 5:     F = A[0]*EZ + A[3] + A[4]*X
  61. 6:     F = A[0]*EZ + A[3] + A[4]*X + A[5]*X^2 ;FUNCTIONS.
  62.     ENDCASE
  63.  
  64.     IF N_PARAMS(0) LE 3 THEN RETURN ;NEED PARTIAL?
  65. ;
  66.     PDER = FLTARR(N_ELEMENTS(X),n) ;YES, MAKE ARRAY.
  67.     PDER[0,0] = EZ        ;COMPUTE PARTIALS
  68.     if a[2] ne 0. then PDER[0,1] = A[0] * EZ * Z/A[2]
  69.     PDER[0,2] = PDER[*,1] * Z
  70.     if n gt 3 then PDER[*,3] = 1.
  71.     if n gt 4 then PDER[0,4] = X
  72.     if n gt 5 then PDER[0,5] = X^2
  73.     RETURN
  74. END
  75.  
  76.  
  77.  
  78. Function Gaussfit, x, y, a, NTERMS=nt, ESTIMATES = est
  79. ;+
  80. ; NAME:
  81. ;    GAUSSFIT
  82. ;
  83. ; PURPOSE:
  84. ;     Fit the equation y=f(x) where:
  85. ;
  86. ;         F(x) = A0*EXP(-z^2/2) + A3 + A4*x + A5*x^2
  87. ;             and
  88. ;        z=(x-A1)/A2
  89. ;
  90. ;    A0 = height of exp, A1 = center of exp, A2 = sigma (the width).
  91. ;    A3 = constant term, A4 = linear term, A5 = quadratic term.
  92. ;    Terms A3, A4, and A5 are optional.
  93. ;     The parameters A0, A1, A2, A3 are estimated and then CURVEFIT is 
  94. ;    called.
  95. ;
  96. ; CATEGORY:
  97. ;    ?? - fitting
  98. ;
  99. ; CALLING SEQUENCE:
  100. ;    Result = GAUSSFIT(X, Y [, A])
  101. ;
  102. ; INPUTS:
  103. ;    X:    The independent variable.  X must be a vector.
  104. ;    Y:    The dependent variable.  Y must have the same number of points
  105. ;        as X.
  106. ; KEYWORD INPUTS:
  107. ; KEYWORD INPUTS:
  108. ;    ESTIMATES = optional starting estimates for the parameters of the 
  109. ;        equation.  Should contain NTERMS (6 if NTERMS is not
  110. ;        provided) elements.
  111. ;    NTERMS = Set NTERMS to 3 to compute the fit: F(x) = A0*EXP(-z^2/2).
  112. ;       Set it to 4 to fit:  F(x) = A0*EXP(-z^2/2) + A3
  113. ;       Set it to 5 to fit:  F(x) = A0*EXP(-z^2/2) + A3 + A4*x
  114. ;
  115. ; OUTPUTS:
  116. ;    The fitted function is returned.
  117. ;
  118. ; OPTIONAL OUTPUT PARAMETERS:
  119. ;    A:    The coefficients of the fit.  A is a three to six
  120. ;        element vector as described under PURPOSE.
  121. ;
  122. ; COMMON BLOCKS:
  123. ;    None.
  124. ;
  125. ; SIDE EFFECTS:
  126. ;    None.
  127. ;
  128. ; RESTRICTIONS:
  129. ;    The peak or minimum of the Gaussian must be the largest
  130. ;    or smallest point in the Y vector.
  131. ;
  132. ; PROCEDURE:
  133. ;    The initial estimates are either calculated by the below procedure
  134. ;    or passed in by the caller.  Then the function CURVEFIT is called
  135. ;    to find the least-square fit of the gaussian to the data.
  136. ;
  137. ;  Initial estimate calculation:
  138. ;    If the (MAX-AVG) of Y is larger than (AVG-MIN) then it is assumed
  139. ;    that the line is an emission line, otherwise it is assumed there
  140. ;    is an absorbtion line.  The estimated center is the MAX or MIN
  141. ;    element.  The height is (MAX-AVG) or (AVG-MIN) respectively.
  142. ;    The width is found by searching out from the extrema until
  143. ;    a point is found less than the 1/e value.
  144. ;
  145. ; MODIFICATION HISTORY:
  146. ;    DMS, RSI, Dec, 1983.
  147. ;    DMS, RSI, Jun, 1995, Added NTERMS keyword.  Result is now float if 
  148. ;                Y is not double.
  149. ;    DMS, RSI, Added ESTIMATES keyword.
  150. ;-
  151. ;
  152. on_error,2                      ;Return to caller if an error occurs
  153. csave = !c
  154. if n_elements(nt) eq 0 then nt = 6
  155. if nt lt 3 or nt gt 6 then $
  156.    message,'NTERMS must have values from 3 to 6.'
  157. n = n_elements(y)        ;# of points.
  158. s = size(y)
  159.  
  160. if n_elements(est) eq 0 then begin    ;Compute estimates?
  161.     c = poly_fit(x,y,1,yf)        ;Fit a straight line
  162.     yd = y - yf
  163.     if s[s[0]+1] ne 5 then begin    ;If Y is not double, use float
  164.     yd = float(yd) & c = float(c) & endif
  165.  
  166.     ymax=max(yd) & xmax=x[!c] & imax=!c    ;x,y and subscript of extrema
  167.     ymin=min(yd) & xmin=x[!c] & imin=!c
  168.  
  169.     if abs(ymax) gt abs(ymin) then i0=imax else i0=imin ;emiss or absorp?
  170.     i0 = i0 > 1 < (n-2)        ;never take edges
  171.     dy=yd[i0]            ;diff between extreme and mean
  172.     del = dy/exp(1.)        ;1/e value
  173.     i=0
  174.     while ((i0+i+1) lt n) and $    ;guess at 1/2 width.
  175.     ((i0-i) gt 0) and $
  176.     (abs(yd[i0+i]) gt abs(del)) and $
  177.     (abs(yd[i0-i]) gt abs(del)) do i=i+1
  178.     a = [yd[i0], x[i0], abs(x[i0]-x[i0+i])]
  179.     if nt gt 3 then a = [a, c[0]]    ;estimates
  180.     if nt gt 4 then a = [a, c[1]]
  181.     if nt gt 5 then a = [a, 0.]
  182. endif else begin
  183.     if nt ne n_elements(est) then message, 'ESTIMATES must have NTERM elements'
  184.     a = est
  185. endelse
  186.  
  187. !c=csave            ;reset cursor for plotting
  188. return,curvefit(x,y,replicate(1.,n),a,sigmaa, $
  189.         function_name = "GAUSS_FUNCT") ;call curvefit
  190. end
  191.